#!/usr/bin/R
T=as.numeric(try(system("ls simulated.txt*weights | wc -l",intern = TRUE))) ### how many generations we picked up
d<-read.table("analytical.txt.f1");
dims=ncol(d) ## how many dimensions of data 
var<-vector(); for (i in (1:ncol(d))){var<-c(var,(d[i+1,i]))} ## true variance of each dimension (mean assumed to be 0)
scale=10; # scale of the plot, this should be fine

system("echo '<html><body>' > visualize.html");
T=T-1;
for (dim in (1:dims)) {

emperical<-read.table("emperical.txt")[,dim];  ## here is the NUMERICAL data
xbar=mean(emperical); 
s=sd(emperical); 
n=length(emperical);  

u=d[1,dim]; t=sqrt(var[dim]);  ## this is the ANALYTICAL data....
mean=u*1/t^2/(1/t^2+n/s^2)+xbar*n/s^2/(1/t^2+n/s^2); 
sd=sqrt(((t^2*s^2/n)/(t^2+s^2/n))^2); 

plotlims=seq(mean-scale*sd,mean+scale*sd,length.out=1000); 

for (i in (0:T)) {
	png(paste("simulated.txt",i,".png",sep='')); ##mean-scale*sd,mean+scale*sd
	hist(read.table(paste("simulated.txt",i,sep=''))[,dim],xlab=paste("u",dim,sep=""),ylab=paste("density of (u",dim,"1|x)",sep=""),xlim=c(floor(mean-scale*sd),ceiling(mean+scale*sd)),ylim=c(0,dnorm(mean,mean,sd)),prob=T,main=paste("p(u",dim,"|x) simulation:",i+1),breaks=1000); 
	lines(plotlims,dnorm(plotlims,mean=mean,sd=sd),col="red"); 
	dev.off();
}; 

system(paste("convert -delay 50 -loop 0 simulated.txt{0..",T,"}.png simulation",dim,".gif",sep='')); ##requires ImageMagick



system(paste("echo '<img src=\"simulation",dim,".gif\">' >> visualize.html",sep=""));
}

for (i in (0:T)) {
        png(paste("simulated.txt",i,"weights.png",sep=''));

	if (i>1) {
        	plot(read.table(paste("simulated.txt",i,"weights",sep=''))[,1],main=paste("weighting parameters, simulation:",i+1),ylim=c(0,2/1000),type="n",sub=ks.test(read.table(paste("simulated.txt",i-1,"weights",sep=''))[,1],read.table(paste("simulated.txt",i,"weights",sep=''))[,1])[[2]]); 
	}
	else {
		plot(read.table(paste("simulated.txt",i,"weights",sep=''))[,1],main=paste("weighting parameters, simulation:",i+1),ylim=c(0,2/1000),type="n" );
	}
	lines(read.table(paste("simulated.txt",i,"weights",sep=''))[,1]);
	lines(sort(read.table(paste("simulated.txt",i,"weights",sep=''))[,1]),col="red");
        dev.off();

        png(paste("simulated.txt",i,"weights-hist.png",sep=''));
        hist(read.table(paste("simulated.txt",i,"weights",sep=''))[,1],main=paste("weighting parameters, simulation:",i+1),ylim=c(0,1000),xlim=c(0,2/1000));
        dev.off();

}
system(paste("convert -delay 50 -loop 0 simulated.txt{0..",T,"}weights.png simulation-weights.gif",sep=''));
system(paste("convert -delay 50 -loop 0 simulated.txt{0..",T,"}weights-hist.png simulation-weights-hist.gif",sep=''));

system(paste("echo '<img src=\"simulation-weights.gif\">' >> visualize.html",sep=""));
system(paste("echo '<img src=\"simulation-weights-hist.gif\">' >> visualize.html",sep=""));



system("echo '</body></html>' >> visualize.html")
system("rm -v simulated.txt*.png");
